---
title: 'Latta et al 2016 Costa Rica: Species Richness Regressions'
author: "brouwern@gmail.com"
date: "December 22, 2016"
output:
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE)
```

This code carries out analyses related to regression modeling of species richness.  The code has not been edited for clarity or to remove apsects that did not appear in the final manuscript.


# Preliminaries

## Load data

```{r, message=FALSE, warning=FALSE}
#setwd
my.dir <-  "C:/Users/lisanjie2/Dropbox/Latta_et_al_2017_PeerJ_FINAL"

setwd(my.dir)

#load data
CR.long <- read.csv("./DATA/Latta_et_al_mistnet_long.csv")

```

## Packages

```{r, message=FALSE, warning=FALSE}
library(reshape2)
library(ggplot2)
library(doBy)
library(lme4)

# for calculating combined coefficients from interaction models
library(multcomp)  

#for display()
library(arm)
```

## Bolkers prediction functions

```{r}
# from https://rpubs.com/bbolker/glmmchapter
source("./functions/Bolker_easyPredCI_12222016.R")
```



# Species Richness: total

## Calcualte species richness

Make pres-abs
```{r}

CR.melt.all.long$pres.abs <- ifelse(CR.melt.all.long$N > 0,1,0)

```

### Overall species richness
Rehsape to calcautle overall species richness
```{r}
# 
spp.rich <- dcast(data = CR.melt.all.long,
                  formula = year + month ~ .,
                  value.var = "pres.abs",
                  fun.aggregate = sum)

names(spp.rich)[3] <- "tot.spp.rich"
```

### Species richness by status

Species richness by status (lat mig vs resident)
```{r}
spp.rich.by.status <- dcast(data = CR.melt.all.long,
                            formula = year + month  + status ~ . ,
                            value.var = "pres.abs",
                            fun.aggregate = sum)
names(spp.rich.by.status)[4] <- "spp.rich.status"
spp.rich.by.status$status <- factor(spp.rich.by.status$status,
                                    levels = c("PR","M"))

```

### Species richnessby habitat (lat mig vs resident)

```{r}
spp.rich.by.hab.prime <- dcast(data = CR.melt.all.long,
                               formula = year + month + hab.prime ~ .,
                               value.var = "pres.abs",
                               fun.aggregate = sum)


names(spp.rich.by.hab.prime)[4] <- "spp.rich"

```


### Species richness by foraging guild2 (lat mig vs resident)

NOTE: Some foraging guilds are rare.  Forage guild2 pools rare groups.
```{r}
spp.rich.by.forage.guilde2 <- dcast(data = CR.melt.all.long,
                                    formula = year + month + forage.guilde2 ~ .,
                                    value.var = "pres.abs",
                                    fun.aggregate = sum)


names(spp.rich.by.forage.guilde2)[4] <- "forage.guilde2.rich"

```

### Species richness by sensitivity to disturbance
```{r}



# 
spp.rich.by.sense <- dcast(data = CR.melt.all.long,
                           formula = year + month + sense2 ~ .,
                           value.var = "pres.abs",
                           fun.aggregate = sum)


names(spp.rich.by.sense)[4] <- "frag.sense.rich"



```

Mean species richness by month
```{r kable}
library(knitr)
sum.out <- summaryBy(tot.spp.rich ~ month, data = spp.rich,
           FUN = c(mean),na.rm =T)
kable(sum.out)
```


## Model overall species richness

### Fit overall species richness models
```{r, message=FALSE, warning=FALSE, error = FALSE}
### OVERALL SPECIES RICHNESS

#set control options to overide erros, which are not playing nice w/rmarkdown for some reason..
cntrl <- lmerControl(check.nobs.vs.nlev = "ignore",check.nobs.vs.rankZ =   
   "ignore",check.nlev.gtreq.5 = "ignore",check.nobs.vs.nRE="ignore",
   check.rankX =    c("ignore"),
   check.scaleX = "ignore",
   check.formula.LHS="ignore",
   check.conv.grad   = .makeCC("warning", tol = 1e-3, relTol = NULL))

#null model
lme.spp.rich.overall.null <- lmer(tot.spp.rich ~  1 + (1|year) , data = spp.rich,control = cntrl)

#month only (intercept)
lme.spp.rich.overall.month <- update(lme.spp.rich.overall.null, .~. + month)

#year
lme.spp.rich.overall.year <- update(lme.spp.rich.overall.null, .~. + year)

#month + year
lme.spp.rich.overall.add <- update(lme.spp.rich.overall.null, .~. +month+ year)

#full model: month*year
lme.spp.rich.overall.full <- update(lme.spp.rich.overall.null, .~. + month*year)

#calcualte means
lme.spp.rich.overall.month.means <- lmer(tot.spp.rich ~  -1+month + (1|year) , data = spp.rich)

display(lme.spp.rich.overall.month.means, detail = T)
```

### Test effects from overall species richness
```{r, message=FALSE, warning=FALSE}
### Month effect
anova(lme.spp.rich.overall.month,
      lme.spp.rich.overall.null)

### year effect
anova(lme.spp.rich.overall.year,
      lme.spp.rich.overall.null)

###
anova(lme.spp.rich.overall.add,
      lme.spp.rich.overall.year)


anova(lme.spp.rich.overall.add,
      lme.spp.rich.overall.year)


### interaction
anova(lme.spp.rich.overall.add,
      lme.spp.rich.overall.full)
```


## Model species group by latitudinal migration behavior

### Fit migrant vs resident species richness models
```{r}
i.Jan <- which(spp.rich.by.status$month == "Jan")
lme.by.status.Jan.null <- lmer(log(spp.rich.status) ~ 1 + (1|year),
                               data = spp.rich.by.status[i.Jan,])
lme.by.status.Jan.yr     <- update(lme.by.status.Jan.null, .~.+scale(year))
lme.by.status.Jan.status <- update(lme.by.status.Jan.null, .~.+status)

lme.by.status.Jan.add  <- update(lme.by.status.Jan.null, .~. + scale(year) + status)
lme.by.status.Jan.full <- update(lme.by.status.Jan.null, .~. + scale(year)*status)

```


### Test migrant vs resident species richness models
```{r, message=FALSE, warning=FALSE}
# year effect
anova(lme.by.status.Jan.yr, 
      lme.by.status.Jan.null)

anova(lme.by.status.Jan.full,
      lme.by.status.Jan.add)
summary(lme.by.status.Jan.full)
```


### Calculate slopes

Use multcomp pacakge to combine coefficient from year + trait*year
```{r}

### Calculate beta, SE and p for scale(year)+scale(year):status
contrast.matrix <- rbind("contrast" = c(0, 1, 0, 1))
mult.comp.working <- glht(lme.by.status.Jan.full       
                          ,linfct =contrast.matrix
                          ,alternative = c("two.sided"))
mult.comp.summary <- summary(mult.comp.working
                             ,test = adjusted(type = "none"))

mult.comp.summary

```



### Plot output of migrant vs. residents model

Set up labels
```{r}
spp.rich.by.status$status2 <- as.character(spp.rich.by.status$status)
spp.rich.by.status$status2 <- ifelse(spp.rich.by.status$status == "PR", "Resident", "Migrant")


```


#### Run model on raw data 
No transformation so that params can be directly interp
```{r}
#model rawdata
lme.by.status.Jan.full.RAWDAT <- lmer(spp.rich.status ~ year*status + (1|year),
                                      data = spp.rich.by.status[i.Jan,],control = cntrl)

display(lme.by.status.Jan.full.RAWDAT)


#predictions from lmer
pred.out <- easyPredCI.lmm(model = lme.by.status.Jan.full.RAWDAT,
                           newdata = lme.by.status.Jan.full.RAWDAT@frame)


lme.by.status.pred.out <- cbind(pred.out, lme.by.status.Jan.full.RAWDAT@frame)

```


#### Plot migrant vs. resident species richenss
Load my ggplot function
```{r}
source("./functions/Fx_ggPlot_sppRichRegxn.R")
```


```{r}


# Plot species richess ~ migration status

gg.spp.rich.mig.status <-gg.spp.rich.regressions.nofacets(gg.dat = lme.by.status.pred.out,
                                 x = "year",
                                 y1 = "predict",
                                 y2 = "spp.rich.status",
                                 color. = "status",
                                 shape. = "status",
                                 linetype. = "status",
                                 fill. = "status",
                                 legend.name = "Migratory Status",
                                 legend.breaks = c("PR",
                                                   "M"),
                                 legend.labels = c("Resident",
                                                   "Migrant"),
                                 alpha. = 0.1) + theme(
  legend.position=c(0.105, .935))
gg.spp.rich.mig.status

```



## Model species richness by habitat

Primary forest vs. secondary forest species
```{r, message=FALSE, warning=FALSE}
# single terms - no time interactions
lme.hab.prime.null  <- lmer(log(spp.rich) ~ 1 +(1|year) + 
                              (1|year:month), 
                            data = spp.rich.by.hab.prime)
lme.hab.prime.month <- lmer(log(spp.rich) ~ month +(1|year) +(1|year:month), data = spp.rich.by.hab.prime)
lme.hab.prime.hab   <- lmer(log(spp.rich) ~ hab.prime +(1|year) +(1|year:month), data = spp.rich.by.hab.prime)
lme.hab.prime.yr    <- lmer(log(spp.rich) ~ scale(year) +(1|year) +(1|year:month), data = spp.rich.by.hab.prime)
lme.hab.prime.month.hab.add    <- lmer(log(spp.rich) ~hab.prime + month +(1|year) +(1|year:month), data = spp.rich.by.hab.prime)
lme.hab.prime.month.hab.mult    <- lmer(log(spp.rich) ~hab.prime*month +(1|year) +(1|year:month), data = spp.rich.by.hab.prime)

### Test for habitat difference
anova(lme.hab.prime.hab,
      lme.hab.prime.null)

anova(lme.hab.prime.month.hab.mult,
      lme.hab.prime.month.hab.add)


#model change in primary forest species over time
lme.hab.prime.yrXhab    <- lmer(log(spp.rich) ~ hab.prime*scale(year) +(1|year) +(1|year:month), data = spp.rich.by.hab.prime)

lme.hab.prime.best.drop1    <- lmer(log(spp.rich) ~ month*hab.prime + hab.prime + scale(year) +(1|year) +(1|year:month), data = spp.rich.by.hab.prime)

lme.hab.prime.BEST    <- lmer(log(spp.rich) ~ 1 + month*hab.prime + 
                                hab.prime*scale(year) +
                                (1|year) +(1|year:month), data = spp.rich.by.hab.prime)



lme.hab.prime.ADD.ADD    <- lmer(log(spp.rich) ~ month*scale(year) +hab.prime*scale(year)+(1|year) +(1|year:month), data = spp.rich.by.hab.prime)


lme.hab.prime.FULL    <- lmer(log(spp.rich) ~ month*hab.prime*scale(year) +(1|year) +(1|year:month), data = spp.rich.by.hab.prime)


### Habitat effect 
anova(lme.hab.prime.hab,lme.hab.prime.null)

```


Best model
```{r}
### BEST MODEL:Year*habitat.prime ###
summary(lme.hab.prime.BEST)

anova(lme.hab.prime.BEST, lme.hab.prime.best.drop1)


#check in CAR

library(car)
Anova(lme.hab.prime.BEST, type = "III")
```


### Calcualte slopes
```{r}
## calcualte slope 
library(multcomp)
summary(lme.hab.prime.BEST)

fixef(lme.hab.prime.BEST)
### Calculate beta, SE and p for scale(year)+ hab.primeS:scale(year)
contrast.matrix <- rbind("contrast" = c(0, 0, 0, 1, 0, 1))
mult.comp.working <- glht(lme.hab.prime.BEST       
                          ,linfct =contrast.matrix
                          ,alternative = c("two.sided"))
mult.comp.summary <- summary(mult.comp.working
                             ,test = adjusted(type = "none"));mult.comp.summary

```


### Plot

#### Refit model on untransformed data
```{r, message=FALSE, warning=FALSE}
lme.hab.prime.BEST.RAWDAT    <- lmer(spp.rich ~ 1 + month*hab.prime + 
                                       hab.prime*year +
                                       (1|year) +(1|year:month), 
                                     data = spp.rich.by.hab.prime,control = cntrl)
```

#### predictions from lmer
```{r}

pred.out <- easyPredCI.lmm(model = lme.hab.prime.BEST.RAWDAT,
                           newdata = lme.hab.prime.BEST.RAWDAT@frame)

lme.hab.prime.pred.out <- cbind(pred.out, lme.hab.prime.BEST.RAWDAT@frame)

lme.hab.prime.pred.out$month <- as.character(lme.hab.prime.pred.out$month)
lme.hab.prime.pred.out$month <- ifelse(lme.hab.prime.pred.out$month == "Aug","August","January")
spp.richness.prim.habitat <- gg.spp.rich.regressions(gg.dat = lme.hab.prime.pred.out,
                                                     color. = "hab.prime",
                                                     shape. = "hab.prime",
                                                     linetype. = "hab.prime",
                                                     fill. = "hab.prime",
                                                     legend.name = "Primary Habitat",
                                                     legend.breaks = c("F","S"),
                                                     legend.labels = c("Mature forest",
                                                                       "Secondary forest"),
                                                     alpha. = 0.1) + theme(
  legend.position=c(0.105, .915))


spp.richness.prim.habitat


```



## Model species richnes by foraging guild

Species Richness by guilde (F vs S)
```{r, message=FALSE, warning=FALSE}

mod.by.guild2.best <- lmer(log(forage.guilde2.rich) ~
                            scale(year)*forage.guilde2 + 
               forage.guilde2*month +
               (1|year) +
               (1|year:month), #family = poisson,
             data = spp.rich.by.forage.guilde2[,])


mod.by.guild2.best.drop1 <- lmer(log(forage.guilde2.rich) ~
                                  scale(year) + forage.guilde2 + 
                                  forage.guilde2*month +
                                  (1|year) +
                                  (1|year:month), #family = poisson,
                                data = spp.rich.by.forage.guilde2[,])


```

Test models
```{r, message=FALSE, warning=FALSE}
anova(mod.by.guild2.best,mod.by.guild2.best.drop1)

display(mod.by.guild2.best)
```


### Plot species richness ~ guild

#### Refit model on raw data
```{r}
lme.forage.guilde.BEST.RAWDAT  <- lmer(forage.guilde2.rich ~
                                         scale(year)*forage.guilde2 + 
                                         forage.guilde2*month +
                                         (1|year) +
                                         (1|year:month), #family = poisson,
                                       data = spp.rich.by.forage.guilde2[,],control = cntrl)

```

Predictons from lmer
```{r}

#predictions from lmer
pred.out <- easyPredCI.lmm(model = lme.forage.guilde.BEST.RAWDAT,
                           newdata = lme.forage.guilde.BEST.RAWDAT@frame)

lme.forage.guilde.pred.out <- cbind(pred.out, lme.forage.guilde.BEST.RAWDAT@frame)

lme.forage.guilde.pred.out$month <- as.character(lme.forage.guilde.pred.out$month)
lme.forage.guilde.pred.out$month <- ifelse(lme.forage.guilde.pred.out$month == "Aug","August","January")

```

Plot
```{r, message=FALSE, warning=FALSE}


### Plot species richness ~ Foraging guild
library(ggplot2)
summary(lme.forage.guilde.pred.out)

lme.forage.guilde.pred.out$forage.guilde2 <- gsub("O","Omnivore",lme.forage.guilde.pred.out$forage.guilde2 )

lme.forage.guilde.pred.out$forage.guilde2 <- gsub("I","Insectivore",lme.forage.guilde.pred.out$forage.guilde2 )

lme.forage.guilde.pred.out$forage.guilde2 <- gsub("F/N","Frugivore/Nectarivore",lme.forage.guilde.pred.out$forage.guilde2 )


gg.spp.rich.guild.habitat <- gg.spp.rich.regressions(gg.dat = lme.forage.guilde.pred.out,
                                              color. = "forage.guilde2",
                                              shape. = "forage.guilde2",
                                              linetype. = "forage.guilde2",
                                              fill. = "forage.guilde2",
                                              legend.name = "Primary\nHabitat",
                                              legend.breaks = c("F/N",
                                                                "I",
                                                                "O"),
                                              legend.labels = c("Frug./Nectivore",
                                                                "Insectivore",
                                                                "Omnivore"),
                                              alpha. = 0.1,
                                              y2 = "forage.guilde2.rich")+
                                              facet_grid(forage.guilde2~month) #change facets


gg.spp.rich.guild.habitat





```



## Species richness by sensitivity to fragmentation

Species Richness by frag/dist sensitivity

```{r}

frag.sense.sense <- lmer(log(frag.sense.rich) ~ scale(year)*sense2 + month + #sense2*month +
               (1|year) +
               (1|year:month), #family = poisson,
             data = spp.rich.by.sense)

frag.sense.best <- lmer(log(frag.sense.rich) ~ scale(year)*sense2 + month + #sense2*month +
               (1|year) +
               (1|year:month), #family = poisson,
             data = spp.rich.by.sense)
frag.sense.drop1 <- lmer(log(frag.sense.rich) ~ scale(year) + sense2 + month + #sense2*month +
                          (1|year) +
                          (1|year:month), #family = poisson,
                        data = spp.rich.by.sense)
anova(frag.sense.best,frag.sense.drop1)

```


### Plot species richness ~ primary habitat

#### Refit model
```{r}

lme.fragsense.BEST.RAWDAT <- lmer(frag.sense.rich  ~ year*sense2 + month + #sense2*month +
                                    (1|year) +
                                    (1|year:month), #family = poisson,
                                  data = spp.rich.by.sense,control = cntrl)

```


### predictions from lmer

```{r}
#
pred.out <- easyPredCI.lmm(model = lme.fragsense.BEST.RAWDAT,
                           newdata = lme.fragsense.BEST.RAWDAT@frame)

lme.fragsense.pred.out <- cbind(pred.out, lme.fragsense.BEST.RAWDAT@frame)

lme.fragsense.pred.out$month <- as.character(lme.fragsense.pred.out$month)
lme.fragsense.pred.out$month <- ifelse(lme.fragsense.pred.out$month == "Aug","August","January")
```


### Plot
```{r}




gg.spp.rich.frag.sense <- gg.spp.rich.regressions(gg.dat = lme.fragsense.pred.out,
                                                  color. = "sense2",
                                                  shape. = "sense2",
                                                  linetype. = "sense2",
                                                  fill. = "sense2",
                                                  legend.name = "Fragmentation Sensitivity",
                                                  legend.breaks = c("H/M","L"),
                                                  legend.labels = c("Sensitive",
                                                                    "Insensitive"),
                                                  alpha. = 0.1,
                                                  y2 = "frag.sense.rich") +  theme(
  legend.position=c(0.1575205, .915)) 


gg.spp.rich.frag.sense 
```










# Save output
```{r}

date.str <- format(Sys.time(), "%a%b%d%Y")

workspace.str <- paste("workspace_all_spp_richens",date.str,
                       ".RData", 
                       sep = "")

save.image(file = workspace.str)
```









# Session info
```{r session_info, include=TRUE, echo=TRUE, results='markup'}
devtools::session_info()
```